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Objectives 

Microspacecraft  are  currently  considered  by  Air  Force,  NASA  and  industry  for  a  variety  of  applications. 
Characteristics  of  a  microspacecraft  include  dimensions  less  than  1 0  cm,  weight  less  than  1  kg,  and  power 
below  10  W.  These  requirements  clearly  necessitate  the  need  for  microtechnologies  and  among  them  a 
critical  one  is  micropropulsion.  The  mathematical  and  computational  work  investigated  new  issues  that 
arise  in  simulating  gas  and  plasma  microthrusters  that  are  an  order  of  magnitude  smaller  than  the  existing 
state-of-the  art.  Flows  in  these  microthrusters  exhibit  distinct  physical  characteristics  that  make 
computational  modeling  challenging.  Current  simulation  technologies  are  based  on  continuum 
mathematical  models  and  assumptions  that  break  down  in  microscales  and,  therefore,  cannot  treat  these 
microflows  comprehensively.  In  addition,  unlike  the  classical  aerodynamic  flows  that  are  governed  by 
predominantly  hyperbolic  behavior,  flows  in  microthrusters  are  strongly  viscous  and,  thus,  exhibit 
parabolic  behavior.  Plasma  microthrusters  introduce  an  additional  complexity  in  their  analysis  due  to  the 
overlap  of  electrodynamic  and  gasdynamic  scales. 

This  research  aimed  overall  at  a  unified  mathematical  formulation  and  numerical  discretization  of 
multi-species,  partially  ionized  plasma  micro-flows  in  non-equilibrium  and  involved  closed  interaction 
between  two  research  groups.  The  group  lead  by  Prof.  Gatsonis  at  WPI  developed  a  seamless  Particle-In- 
Cell/Monte  Carlo  methodology  and  the  group  lead  by  Prof.  Kamiadakis  at  Brown  University  developed  a 
semi-continuous  spectral/hp  methodology  based  on  the  multi-species,  compressible,  viscous-MHD 
equations.  The  developed  methodologies  were  applied  to  the  modeling  of  the  micro-Pulsed  Plasma 
Thruster,  other  micropropulsion  thrusters,  and  related  devices. 

The  major  accomplishments  are  outlined  bellow: 

•  Developed  an  adaptive  unstructured  grid  generator  with  an  enhanced  optimization  capability 
(Kovalev,  2000;  Hammel,  2002) 

•  Developed  an  efficient  particle-transfer  algorithm  during  grid  adaptation. 

•  Implemented  neutral  and  charged  particle  transport  algorithms  on  unstructured  grids. 

•  Implemented  elastic  neutral-neutral,  elastic  ion-neutral  and  charge  exchange  ion-neutral  collisions  via 
the  no-time-counter  method. 

•  Implemented  the  Larsen-Borgnakee  model  for  energy  redistribution  for  rotational  and  vibrational 
degrees  of  freedom. 

•  Implemented  a  finite-volume  method  for  Poisson’s  equation  that  utilized  the  Voronoi-Delauney  mesh 
and  implemented  a  GMRES  solver  to  the  discrete  Poisson  equation. 

•  Performed  simulations  of  non-ionized  flows  in  micronozzles  as  well  as  MEMS-based  micronozzles 
(Hammel  et  al.,  2001,  Hammel,  2002). 

•  Performed  simulations  of  plasma  flows  in  electric  micropropulsion  devices  such  as  the  micro  Field 
Emission  Array  (Gatsonis  and  Spirkin,  2002)  and  the  micro  Retarding  Potential  Analyzer  (Spirkin 
and  Gatsonis,  2003). 

•  Investigated  high-order  particle/force  weighting  and  evaluated  of  fluctuations  in  the  unstructured 
simulation  of  collisionless,  fully-  ionized  plasmas  (Spirkin  and  Gatsonis,  2003). 

•  Assisted  Brown  University  in  the  development  of  an  ablation  model  for  micro-PPT  simulations. 
Simulations  also  are  underway  that  couple  the  continuous  results  from  the  multi-fluid  code  (Brown 
University)  with  the  DSMC/PIC  plume  simulations. 

Accomplishments/New  Findings 
1.  Unstructured  3D  PIC/DSMC  Methodology 
1.1  Grid  Generation  and  Adaptation 

An  unstructured  grid  generator  has  been  developed  that  provides  three-dimensional  meshes  for  arbitrary 
geometric  configurations  and  adapts  existing  meshes  according  to  the  preliminary  solution  obtained  on  an 
initial  grid  (Kovalev,  2000;  Hammel,  2002).  The  unstructured  mesh  generator  is  based  on  Watson’s 
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(1981)  incremental  node  insertion  method,  which  uses  properties 
of  the  Delaunay  triangulation.  An  initial  mesh  is  required  for 
Watson’s  method,  in  order  to  have  a  domain  where  point 
insertion  to  begin.  The  initial  mesh  chosen  is  a  cube  divided  into 
six  tetrahedra.  After  the  initial  mesh  is  generated,  the  source 
geometry  is  inserted  into  the  domain.  The  underlying  sizing 
function  -  defined  by  the  source  geometry  -  requires  the  interior 
of  the  grid  to  be  enriched  with  nodes  to  the  specified  density.  For 
this  purpose  the  algorithm  by  Borouchaki  and  George  (1997)  was 
implemented  and  extended  to  three  dimensions.  In  this 
algorithm,  the  characteristic  distance  between  nodes  is  specified 
for  each  grid  node  based  on  mean-free  path  or  the  Debye  length. 
Every  existing  edge  of  the  grid  is  divided  into  a  number  of  new 
prospective  nodes,  so  that  the  new  resulting  edge  segments  vary 
in  length  gradually  between  the  h-values  of  the  edge  vertices.  The  prospective  nodes  are  filtered  in  order 
to  satisfy  the  spacing  and  grid  quality  criteria.  Nodes  falling  too  close  to  existing  nodes  are  eliminated. 
Nodes  that  worsen  grid  quality  as  specified  by  the  lowest  dihedral  angle  in  a  set  of  cells  are  also 
discarded.  The  nodes  that  are  not  rejected  are  inserted  into  the  grid  via  Watson’s  algorithm.  The  edge 

Minimum  Dihedral  Angle  Distribution 


Figure  1.  Surface  discretization  of 
a  micronozzle. 


Figure  2.  Cross-section  of  grid  used  in  the  simulation  of  a  plenum/nozzle/plume  flow  of  a  microthruster. 
(Left)  Details  of  the  plenum/nozzle  grid  (Center)  Overall  simulation  area  (Right)  Dihedral  angle 
distribution  based  on  the  heuristic  optimization  criterion.  (Hammel  et  al.,  2001) 

division  process  is  repeated  while  new  nodes  are  inserted.  An  example  of  the  surface  discretization  is 
shown  in  Figure  1  and  of  the  grid-generator  in  Figure  2.  A  heuristic  quality  improvement  procedure  has 
been  also  implemented  based  on  a  user-defined  minimum  dihedral  angle.  Figure  2  shows  the  improved 
distribution  of  dihedral  angles  with  the  application  of  optimization. 

Adaptation  is  necessary  for  quality  PIC,  DSMC  or  PIC/DSMC  computations  in  absence  of  a  priori 
knowledge  of  the  flow  because  the  cell  size  must  be  of  the  order  of  the  Debye  length  or  the  mean  free 
path,  which  is  an  unknown  local  flow  characteristic.  Also,  regions  in  which  high  gradients  of  flow 
variables  occur  have  to  be  refined  in  order  to  provide  a  quality  solution,  and  regions  with  very  small 
gradients  are  coarsened  to  save  computational  time.  The  developed  grid  generator  and  the  flow  solver 
interact  with  each  other,  so  that  the  adaptation  is  performed  automatically  after  the  initial  solution  is 
obtained.  In  addition  we  developed  algorithms  for  the  fast-transfer  of  data  between  adaptation. 


1.2  Particle  Loading,  Injection  and  Motion 

The  general  procedures  for  loading  and  injection  used  in  this  work  follows  Birdsall  and  Langdon  (1991) 
and  Bird  (1994).  The  equation  that  describes  the  position  of  a  neutral  particle  is  updated  from  the  post¬ 
collision  velocity  according  to 
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dr 

dt 

The  equations  that  describe  the  motion  of  a  charged  particle  are 


dr 


dv 


■  =  v,  —  =  —  (E  +  vxB) 
m 


(2) 


dt  dt  ric  t 

Integration  is  based  on  a  leapfrog  scheme  following  a  method  developed  by  Boris,  (1970).  Particles  are 

moved  between  adjacent  tetrahedral 
cells  using  a  particle-tracing 
technique.  Following  Figure  3,  the 
intersection  of  a  particle  with  initial 
position  r0  and  velocity  v  with  face 
ABC  is  ^ 

r„  +  vr  =  aAB  -I-  /3A  C  (3) 
where,  x  is  the  time  elapsed  in 
moving  from  the  initial  point  to  the 
plane  defined  by  points  A,  B,  and  C; 
a  and  p  define  the  point  of 
intersection  in  the  skewed  coordinate 
system  of  face  ABC.  If  the  time  x  is 
negative,  an  intersection  does  not 
occur.  If  a  or  p  are  less  than  zero  or 
greater  than  unity,  an  intersection  does  not  occur  within  the  face.  If  the  sum  of  a  and  p  is  greater  than 
unity,  then  the  intersection  occurs  outside  the  face.  Particle-tracking  algorithms  are  particularly  slow  in 
unstructured  grids  and  fast  particle-search  algorithms  are  currently  under  investigation. 


Figure  3.  Tethrahedral  cell. 


Figure  4.  Delaunay- 
Voronoi  mesh  in  2D. 


1.3  Collision  Methodology 

The  particle  motion  and  collisions  are  decoupled  in  our  implementation,  following  the  standard 
DSMC  procedure.  Elastic  collisions  are  based  on  the  no  time  counter  (NTC)  methodology  (Bird,  1994) 
with  molecular  cross-section  given  by  the  variable  hard  sphere  (VHS)  model.  Simulation  of  energy 
transfer  to  and  from  rotational  degrees  of  freedom  is  based  on  the  phenomenological  model  of  Larsen- 
Borgnakke.  Vibrational  degrees  of  freedom  are  modeled  as  discrete,  evenly  spaced  energy  levels  using 
the  harmonic  oscillator  model.  Vibrational  energy  redistribution  is  performed  prior  to  rotational  energy 
redistribution  in  the  collision  process,  following  Bird.  If  vibrational  energies  become  high  enough,  the 
spacing  between  levels  lessens  and  an  anharmonic  model  must  be  used  which  may  include  dissociation  as 
outlined  by  Bird  (Hammel,  2002). 


1.4  Poisson’s  Equation  Solver 

Solution  of  Poisson’s  equation 

V2$  =  -£ qtnt  +  qene  eo  (4) 

t=l  / 

is  based  on  a  finite  volume  formulation  that  takes  advantage  of  the  Voronoi  dual  of  the  Delaunay 
triangulation.  The  Voronoi  cell  corresponding  to  each  Delaunay  node  contains  the  set  off  points  closer  to 
the  point  that  any  other,  the  facets  of  the  Voronoi  cell  are  orthogonal  to  the  lines  joining  the  tetrahedral 
nodes  as  shown  in  Figure  4.  For  a  node-centered  finite-volume  scheme  with  finite  volume  i  associated 
with  node  i  with  a  number  of  faces  Np  the  semi-discrete  equation  form  of  Gauss’s  law  is 


k= i  >,*  l  v 


dn 


Qi 


(5) 


i,Jb 


In  the  above,  Q,  is  the  total  charge  enclosed  by  volume  i,  Ai  k  is  the  area  of  the  face  k  of  associated  with 
volume  i,  and  the  summation  is  over  all  the  faces  of  the  finite  volume  as  shown  in  Figure  4  for  a  2-D  case. 
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Also,  Ai,it=(A  n  )i(k  is  the  magnitude  of  the  area  of  the  face  k  for  the  volume  associated  with  node  i 
multiplied  by  the  unit  outward  normal  vector.  An  approximation  of  the  electric  flux  into  a  cell 
containing  node  i  across  a  Voronoi  face  corresponding  to  an  edge  with  nodes  k  and  i  is: 

£v$.(hA)it=^($t-^)  (6) 

k=i  '  A,j t 


In  the  above,  L{  k  is  the  distance  from  node  i  to  node  k.  Summing  over  all  faces  k  of  a  Voronoi  cell 

corresponding  to  node  i,  a  system  of  linear  equations  may  be  formed  assuming  the  charge  inside  the 
volume,  Q{ ,  is  known.  The  method  reduces  to  the  standard  2nd  order  finite-difference  method  on 

Cartesian  meshes.  The  matrix  form  of  this  equation  is: 


Xi 

*1,2 

*1,3 

*  Run 

'V 

<2, 
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*2,3 
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*3,3 

“  *3,  N 

*3 

_  1 
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*;v,tf 
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(7) 


N  is  the  number  of  nodes  in  the  mesh.  The  coefficients  are  determined  by 

N F,i 

Rij  =  for  i  =  j,  Ri .  = - —  if  j  is  adjacent  to  i,  R.  .  =  0  otherwise. 

fct'  Li  k  '  Litj 


Al 

hi 


is  the  ratio  of  the  area  of  the  Voronoi  face  between 


nodes  i  and  j  to  the  distance  between  nodes  i  and  j  if  the 
nodes.  The  system  of  equations  is  solved  using  GMRES. 

In  a  bounded  domain,  a  solution  of  the  Poisson’s 
equation  is  specified  by  piece-wise  Dirichlet  and 
Neumann  boundary  conditions.  For  the  boundary  shown 
on  the  Figure  5,  Node  1  is  a  Dirichlet  boundary  with 
potential  <J>0 ,  Node  2  is  a  node  on  a  Neumann  boundary 


with  inward  flux  EN  2AN  2 .  Since  the  boundaries  of  the 

Delaunay  mesh  are  forced  to  coincide  with  the  boundaries  of  the  mesh,  the  implementation  of  boundary 
conditions  is  straightforward.  In  the  case  of  Dirichlet  boundary  condition  the  voltage  is  placed  on  the 
right  side  hand  of  the  matrix  and  the  corresponding  row  zeroed  with  a  one  placed  on  the  diagonal.  Fluxes 
in  the  Neumann  boundary  conditions  are  added  to  the  flux  formulation  for  the  Voronoi  cell  corresponding 
to  the  boundary  node,  the  value  of  the  inward  normal  electric  field  multiplied  by  the  boundary  area  is 
added  to  the  right  hand  side  of  the  node  of  interest. 


1.5  Particle  and  Force  Weighting/Interpolation 

Two  weighting  calculations  are  the  part  of  a  PIC  timestep  loop:  accumulation  of  charge  at  grid  nodes 
(weighting)  and  distribution  of  electric  field  to  particle  position  (interpolation).  Numerically  these 
interpolations  are  implemented  in  terms  of  cloud  shape  function  S  or  assignment  function  W  as  in 
Hockney  and  Eastwood  (1999)  or  just  in  terms  of  assignment  function  W ,  which  is  referred  to  as 
weighting  function  (Birdsall  and  Langdon,  1991).  Assignment  functions  of  zero  (NGP)  and  first  order 
(Lagrange  polynomials)  are  used  in  our  PIC/DSMC  code.  Implementation  of  higher  order  schemes  can 
improve  the  quality  of  PIC  simulations.  However,  construction  of  such  schemes  on  3-D  unstructured 
grids  leads  to  complex  analytical  formulations  that  increase  the  computational  time  drastically.  In  order 
to  eliminate  the  self-force  and  conserve  momentum  the  ‘momentum  conservation’  implementation  is 
pursued  where  the  weighting  and  interpolation  function  are  identical. 
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1.6  Numerical  Heating 

An  important  quantity  in  plasma  simulations,  is  the  degree  of  numerical  heating  associated  with  a  non¬ 
physical  increase  in  the  total  energy  of  the  system.  The  numerical  heating  is  described  in  terms  of  the 


Figure  6.  Heating  time  as  a  function  of  At  and  Ar  Figure  7.  Heating  time  as  a  function  of  At  and  Ar 

for  NGP  assignment  scheme.  for  Linear  assignment  scheme 

heating  time  th  ,  which  measures  the  lack  of  energy  conservation  in  the  computational  model  and  is 
infinite  in  real  plasmas.  In  our  analysis  of  the  heating  time  we  use  the  definition  introduced  by  Hockney 
(1971)  as  the  time  that  takes  the  energy  of  system  to  double.  Results  of  the  heating  time  calculations  in 
Figures  6-7  show  that  for  a  specified  \D  /  Ar ,  the  timestep  should  be  chosen  from  the  specific  range  of 
values  to  produce  the  minimum  numerical  error.  From  Figures  6  and  7  it  is  evident  that  the  heating  time 
obtained  when  using  linear  weighting  is  about  10  times  larger  than  the  one  for  the  NGP  weighting.  The 
energy-conserving  implementation  was  also  used  to  calculate  r H ,  but  it  did  not  provide  any  advantage 
over  the  momentum-conserving  scheme.  (Spirkin  and  Gatsonis,  2003) 


Figure  8.  Mach  numbers  in  a  parabolic  15.4-micron 
micronozzle.  Solid  lines  are  from  2D  DSMC 
simulations  of  Piekos  and  Breuer  (1996). 


2.  Gaseous  Micropropulsion  Applications  - 
DSMC 

The  unstructured  DSMC  modules  of  the  code 
were  validated  with  simulation  of  gaseous 
micronozzles  and  comparisons  with  previous 
experimental  and  numerical  results  (Hammel  et  al, 
2001).  Rothe’s  5-mm  diameter  micronozzle 
operating  at  80  Pa  was  simulated  and  results 
compared  favorably  with  the  experiments.  The 
code  was  also  applied  to  the  simulation  of  a 
parabolic  planar  micronozzle  with  a  15.4-micron 
throat.  Results  shown  in  Figure  8  were  compared 
favorably  with  previous  2D  Monte  Carlo 
simulations.  The  Gravity  Probe-B  micronozzle 
was  simulated  in  a  domain  that  includes  the 
injection  chamber  and  plume  region.  Stagnation 
conditions  included  a  pressure  of  7  Pa  and  mass 
flow  rate  of  0.012  mg/s.  The  simulation  shown  in 


Figure  9,  examined  the  role  of  injection  conditions  in  micronozzle  simulations  and  results  are  compared 
with  previous  Monte  Carlo  simulations.  Finally,  the  code  was  applied  to  the  simulation  of  a  34-micron 
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Figure  9.  3D  Unstructured  DSMC  simulation  of  a  combined  plenum/nozzle/plume  flow  of 
a  cold-gas  microthruster.  (Left)Mach  numbers  (Right)  Number  density  In  the  plume 
region  (Hamell  et  al.,  2001). 

throat  MEMS  fabricated  micronozzle.  This  micronozzle  is  planar  in  profile  with  sidewalls  binding  the 
upper  and  lower  surfaces.  The  stagnation  pressure  was  set  at  3.447  kPa  and  represents  an  order  of 
magnitude  lower  pressure  than  used  in  previous  experiments.  The  simulation  demonstrated  the  formation 
of  large  viscous  boundary  layers  in  the  sidewalls.  (Hammel,  2002) 


3.  Electric  Micropropulsion  Applications  -PIC 

In  Gatsonis  and  Spirkin  (2002)  we  presented  the  PIC  implementation  of  our  PIC/DSMC  code  and 
simulations  of  field  emission  array  (FEA)  cathodes  that  can  be  an  integral  part  of  micropropulsion  devices 
(Marrese  et  al.,  2000a).  Electron  emission  in  a  background  plasma  leads  under  certain  conditions  in  the 
formation  of  a  virtual-cathode  that  could  limit  the  operation  of  such  devices.  Our  simulations  expanded 

the  work  of  Marrese  et  al.  (2000b)  by  including 
emission  with  cold  (0.1  eV)  and  hot  (5  eV) 
electrons  into  a  background  5-eV  plasma.  The 
emitted  electron  beam/background  plasma 

coupling  processes  considered  in  our  work  are 

important  not  only  in  the  operation  of  FEA 
cathodes  but  also  in  plume  expansion  as  they 
elucidate  basic  plasma  physics  phenomena. 
Results  in  Figure  10  show  the  formation  of 
virtual  cathode  under  certain  emission 

conditions.  The  simulations  served  also  as 
means  of  validation  of  the  newly  developed  PIC- 
modules  of  our  PIC/DSMC  code. 

In  Spirkin  and  Gatsonis  (2003)  the 

unstructured  PIC  coed  was  used  to  simulate  a 
directional-RPA  currently  under  design  by  the  PI 
to  operate  in  a  high-density  plasma  environment 
(Partridge  et  al.,  2003).  The  simulations  offered 
insights  on  the  plasma  flow  inside  the  device  and 
are  assisting  the  design  process.  The  3D  simulations  results  showed  the  potential  distribution  in  the 
segmented,  micron-scale  channel  as  well  as  the  ion  and  electron  flow  characteristics.  In  an  effort  to 


Case  3 


Case  4 


Figure  10.  Potential  in  front  of  a  field  emission 
array  from  unstructured  3D  PIC  simulations 
(Gatsonis  and  Spirkin,  2002) 
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simulate  the  operation  of  the  sensor,  the  current  on  the  plate  was  evaluated  and  plotted  as  a  function  of  the 
ion  retarding  potential.  The  simulations  presented  the  energy  distribution  of  the  collected  ions. 

In  Wheelock  et  al.  (2003)  the  unstructured  PIC  code  was  implemented  in  the  simulation  of  ion  beam 
neutralization  phenomena.  This  work  is  part  of  collaborative  effort  between  the  PI  and  Dr.  Cooke  of 
AFRL/VSBX,  (Hanscom  AFB).  These  simulations  show  the  dependence  of  the  beam  neutralization  on 
beam  energy  and  neutralization  current.  The  simulations  served  also  as  means  of  validation  of  the  PIC- 
modules  of  our  PIC/DSMC  code  by  comparisons  with  2D  PIC  simulations. 

4.  Development  of  Teflon  Ablation  Model 

This  effort  was  in  collaboration  between  the  PI  and  Professor  Kamiadakis  (Brown  University)  and  led  in 
the  development  of  an  ablation  model  that  is  used  in  the  continuum  simulation  of  PPTs.  The  team 
included  Visiting  International  Scholar  Prof.  Shurzhikov  who  was  hosted  by  WPI  and  Brown  University 
during  the  Summer  of  2001 . 

5.  Axisymmetric  Hybrid  Code 

An  associated  effort  involved  a  hybrid  (particle-fluid)  code  that  was  developed  by  the  PI  with  partial 
support  from  NASA  Glenn  Research  Center  (Gatsonis  and  Yin,  1999).  The  approach  treats  ions  and 
neutrals  as  particles  and  electrons  as  fluid.  We  also  implemented  a  new  conservative  method  for 
weighting  in  axisymmetric  coordinates.  Added  the  full  energy  equation  for  electrons  that  include 
convective,  conductive  and  collisional  terms.  Simulations  of  pulsed  plasma  plumes  were  compared  with 
our  triple  and  quadruple  Langmuir  probe  data  take  in  plumes  of  pulsed  plasma  thrusters.  These  hybrid 
methodologies  validated  in  axisymmetric  coordinates  will  be  transferred  to  the  3D  Unstructured 
DSMC/PIC  code. 
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